package scu.maqiang.tiangong;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;
import java.util.function.Consumer;

public class CTE_T100_QiPing1 {
    public static void main(String[] args) {
        String obj = "DaQiPing1";
        String dir = "C:\\Users\\Administrator\\OneDrive\\TianGongComputation\\T100Computations\\AllMeshes\\Computation_"+ obj + "\\";
        String FluxQP = "FluxData_" + obj + "_";
        String PresQP = "PresData_" + obj + "_";
        String prefixFlux = dir + FluxQP;
        String[] fluxDataFile = {prefixFlux + "100000.txt", prefixFlux + "97500.txt", prefixFlux + "95000.txt", prefixFlux + "92500.txt",
                prefixFlux + "90000.txt", prefixFlux + "87500.txt", prefixFlux + "85000.txt", prefixFlux + "82500.txt",
                prefixFlux + "80000.txt", prefixFlux + "77500.txt" , prefixFlux + "75000.txt"};
        String prefixPres = dir + PresQP;
        String[] presDataFile = {prefixPres + "100000.txt", prefixPres + "97500.txt", prefixPres + "95000.txt", prefixPres + "92500.txt",
                prefixPres + "90000.txt", prefixPres + "87500.txt", prefixPres + "85000.txt", prefixPres + "82500.txt",
                prefixPres + "80000.txt", prefixPres + "77500.txt" , prefixPres + "75000.txt"};
        Mesh3T4 mesh2C = new Mesh3T4().readMesh(dir + "TG" + obj+ "Mesh.mesh", MeshFormat.mesh);
        mesh2C.toTecplot(dir + "My" + obj + "Mesh.dat");
        int dataCount = fluxDataFile.length;
        int nodesInMesh = mesh2C.getNv();
        double[][] hfDataInMesh = new double[dataCount][];
        double[][] prDataInMesh = new double[dataCount][];
        for(int i = 0; i < dataCount; i++) {
            hfDataInMesh[i] = MVO.readArray(fluxDataFile[i]);
            prDataInMesh[i] = MVO.readArray(presDataFile[i]);
        }
        double time = 0.0;
        double dt = 1.0;
        double[] timePoint = {0.0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500};
        double endTime = timePoint[timePoint.length - 1];
        double[] previousFlux = new double[nodesInMesh];
        double[] previousPres = new double[nodesInMesh];
        double[] currentFlux = new double[nodesInMesh];
        double[] currentPres = new double[nodesInMesh];
        Consumer<Double> setCurrentFlux = currenttime -> {
            if(Math.abs(currenttime) < 1.0e-6) {
                Arrays.setAll(currentFlux, i -> hfDataInMesh[0][i]);
                return;
            }
            int location = 0;
            for(int i = 0, n = timePoint.length; i < n; i++) {
                if (currenttime > timePoint[i]) {
                    location++;
                }
            }
            int currentLocation = location;
            Arrays.setAll(currentFlux, i -> {
                double timep = (currenttime - timePoint[currentLocation - 1]) / (timePoint[currentLocation] - timePoint[currentLocation - 1]);
                return (hfDataInMesh[currentLocation][i] - hfDataInMesh[currentLocation - 1][i]) * timep;
            });
        };
        Consumer<Double> setCurrentPressure = currenttime -> {
            if(Math.abs(currenttime) < 1.0e-6) {
                Arrays.setAll(currentPres, i -> prDataInMesh[0][i]);
                return;
            }
            int location = 0;
            for(int i = 0, n = timePoint.length; i < n; i++) {
                if (currenttime > timePoint[i]) {
                    location++;
                }
            }
            int currentLocation = location;
            Arrays.setAll(currentPres, i -> {
                double timep = (currenttime - timePoint[currentLocation - 1]) / (timePoint[currentLocation] - timePoint[currentLocation - 1]);
                return (prDataInMesh[currentLocation][i] - prDataInMesh[currentLocation - 1][i]) * timep;
            });
        };


        String outputHeatFile    = dir + obj + "_heatFile.dat";
        String outputElasticFile = dir + obj + "_ElasticFile.dat";
        String outputStrainFile  = dir + obj + "_StrainFile.dat";
        String outputStressFile  = dir + obj + "_StressFile.dat";

        FES3T41 fsHeat = new FES3T41(mesh2C);
        FES3T42 fsElastic = new FES3T42(mesh2C);
        double kxx = 237, density = 2.7e3, specificHeat = 896.0;
        double elasticE = 6.9e10, elasticNu = 0.28e0, expansion = 2.36e-5;
        double beta = elasticE / (1 - 2 * elasticNu) * expansion;
        double theta = 0.5;
        SRMatrix LMatrix = new SRMatrix(fsHeat.GetNdof());
        SRMatrix RMatrix = new SRMatrix(fsHeat.getNdof());
        fsHeat.assembleHeatThetaScheme(new double[]{kxx, density * specificHeat, dt, theta}, new BVPType[]{BVPType.COMMON, BVPType.CONSIST_MASS}, LMatrix, RMatrix);
        DirectSSolver heatSolver = new DirectSSolver(LMatrix);
        heatSolver.Factorize_Cholesky(1);

        double[] x = new double[fsHeat.getNdof()];
        double[] RHS = new double[fsHeat.getNdof()];
        double[] realTemp = x.clone();
        Arrays.setAll(realTemp,i -> x[i] + 300);
        setCurrentFlux.accept(time);
        Arrays.setAll(previousFlux, i->currentFlux[i]);
        double[] heat_loadFlux = new double[fsHeat.getNdof()];
        fsHeat.assembleFlux(new double[][]{currentFlux}, BVPType.COMMON, heat_loadFlux, 0);
        double[] heat_loadFlux0 = heat_loadFlux.clone();
        mesh2C.toTecplot(outputHeatFile, time, realTemp);
        double[] realHeatFlux = new double[fsHeat.getNdof()];
        while(time < endTime) {
            time += dt;
            setCurrentFlux.accept(time);
            System.out.println("Time: " + time);
            fsHeat.assembleFlux(new double[][]{currentFlux}, BVPType.COMMON, heat_loadFlux, 0);
            double[] Rx = RMatrix.mul(x);
            Arrays.setAll(RHS, i-> Rx[i] + 0.5 * dt * (heat_loadFlux[i] + heat_loadFlux0[i]));
            heatSolver.Solve_Cholesky(RHS, x, 1);
            Arrays.setAll(heat_loadFlux0, i -> heat_loadFlux[i]);
            Arrays.setAll(realTemp, i -> x[i] + 300);
            mesh2C.toTecplot(outputHeatFile, time, realTemp);
        }


    }
}
